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SUMMARY 

We have studied the partition of shear, compressional and kinetic energies in the 
coda of ten earthquakes recorded on a dense array, located at Pinyon Flats Obser- 
vatory (PFO), California. Deformation energies are estimated by measuring finite 
differences of the wavefield components. We have thoroughly studied the validity and 
stability of this technique for the PFO data and obtained reliable measurements in 
the 5-7 Hz frequency band. We observe a clear stabilization of the shear to compres- 
sional (W s /W p ) energy ratio in the coda, with an average value of about 2.8, which 
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is smaller by a factor 2.5 than the expected ratio at the surface of a homogeneous 
Poisson half-space. The ratio between the vertical and horizontal kinetic energies 
(V 2 / H 2 ) can be measured from 5 to 25Hz and shows an abrupt transition from 0.1 
in the 5-10Hz band, to about 0.8 in the 15-25Hz band. These measured values are 
in sharp contrast with the theoretical prediction, around 0.56, for equipartitioned 
elastic waves in a homogeneous half-space. To explain these observations, we have 
developed a theory of equipartition in a layered elastic half-space. Using a rigorous 
spectral decomposition of the elastic wave equation, we define equipartition as a 
white noise distributed over the complete set of eigenfunctions. This definition is 
shown to agree with the standard physical concepts in canonical cases. The theory 
predicts that close to the resonance frequency of a low-velocity layer, the ratio be- 
tween shear and compressional energies strongly decreases. Using a detailed model of 
the subsurface at PFO, this counterintuitive result is found to be in good qualitative 
and quantitative agreement with the observations. Near the resonance frequency of 
the low- velocity structure, the drop of the energy ratios W s /W p and V 2 /H 2 is con- 
trolled by the change of ellipticity of the Rayleigh wave and the large contribution of 
the fundamental Love mode. At higher frequencies, the interplay between Rayleigh 
and Love modes trapped in shallow low-velocity layers is responsible for the abrupt 
increase of the kinetic energy ratio. Our study demonstrates that the partition of 
energy in the seismic coda contains information on the local geological structure. 

Key words: Coda waves, equipartition, layered media, site effect, eigenfunction 
expansion 



1 INTRODUCTION 



Almost forty years ago, lAkil ()1969l ) gave the first explanation of the long tail of seismograms. 
The basic idea is that the diffraction of seismic waves by the heterog eneous crust is at the 
origin of the late arrivals also known as coda waves. Six years later, lAki &: Choueti (119751 ) 
proposed the first quantitative explanation for the algebro-exponential decay of the coda with 
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time. Two models were proposed: (1) The single scattering approximation in which seismic 
waves are deflected at a single scattering site during their propagation between source and 
receiver; (2) The diffusion approximation which basically describes the random walk of the 
seismic energy in the scattering medium. Although conceptually important, these two models 
were not sufficiently accurate to model successfully the s pace time distribution of the energy 



in the coda. This difficulty was partly overcome when 



Wul (|1985l ) introduced the acoustic 



Radiative Transfer Equation (RTE). Mod e conversion s be tween P an d S waves were later 



introduced on a phenomenological basis by 



Zend diaaa) andlSatd p994 ). The RTE takes into 



account arbitrary orders of scattering in a rigorous manner, and describes the behavior of 
intensity in a complex medium beyond the models of single scattering and diffusion. Recently 
the accuracy of th e RTE has been d emonstrated b y comparisons with fu ll-wave simulations 
in both acoustic (jWegler et al.l 120061 ) and elastic (jPrzvbilla et al.l 120061 ) media. While the 
introduction of the RTE in seismology was an important progress, it did not resolve the 



debate initiated by 



Aki fc Chouet 



(|1975l ) on the origin of the coda: is it dominated by single 
or multiple scattering? This question is of both fundamental and practical interest: it is only 
when the propagation regime is known that it becomes possible to separate scattering from 
absorption. 

In the nineties, the RTE f or polar i zed e lastic waves in an i nhom ogeneous medium has 



been derived independently by 



Weaver 


(1990 


) and 


Rvzhik et al. 


(1996) 



(1996). A particularly inter- 



esting outcome of the RTE is the predicted stabilization of the shear (W s ) to compressional 
(W p ) energy ratio (W s /W p ) at large lapse time in the co da, independent of the propertie s 



of the scatterers. This property is known as e quipartition (jWeaver 
and was originally proposed by 



1990 



Rvzhik et al 



19M) 



Weaverl dl982) to describe diffuse waves, based on statisti- 



cal physics arguments. Equipartition implies that while the compressional and shear energy 
decrease exponentially w ith time, their r atio W s /W p tends to a constant. This result was 



numerically confirmed by 



Margerin et al 



(2000) who found that this behavior can occur after 



a few scattering even t s only . A stabilization of the energy ratio W s /W p was observed in Mex- 



ico by 



Shapiro et al.l (|2000i ) and was found to be in excellent quantitative a g reeme nt with an 

(|200lh . The very 



Hennino et al 



equipartition theory for Rayleigh, Love and body waves by 
rapid stabilization of energy ratios observed in Mexico is a clear proof that the propagation of 
short period elastic waves in tectonically active regions is dominated by multiple scattering. 

In this paper, we report the stabilization of the ratio W s /W p , and the ratio between ver- 
tical and horizontal kinetic energies {V 2 /H 2 ), in the coda of short period seismic waves. The 
data were collected during a temporary experiment conducted at the Pinyon Flats Observa- 
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tory (PFO) in California (|Fletcher et al. 



1990 



Vernon et al. 



1998), where a small aperture 



array was installed during six months in 1990. The observed stabilization ratio is around 



2.8 which is about 2.5 times sma 



Poisson half-space (jHennino et al 



ler t han the value expected at the surface of an elastic 
2001). To explain this observation we have developed an 



equipartition theory for elastic waves in a stratified medium. Our approach is based on a spec- 
tral decomposition of the elastic operator developed in section 2, which involves generalized 



eig enfunctions. T 
by iHennino et al 



l is tech nique circumvents the artificial "locked-mode" approximation used 
( 200ll ) and is in much closer agreement with the situation in the Earth. 



ly as a white n oise in the modal space. U sing 
Weaverl (|1990l ) and IHennino et al.1 (|200ll ) are 



Equipartition can then be defined mathematica 
this formal definition, the theoretical results of 
recovered and generalized to stratified media (section 3). The method avoids mode counting 
arguments and treats on the same footing surface and body waves. The data analysis and 
the measurements are presented in section 4. Particular attention has been paid to the stable 



measurements of partial derivatives of the wavefield. In section 5, we show that e c 
theory in a stratified model of the crust at PFO as proposed by 



Fletcher et al 



uipar tition 



(|l99d ) is in 



good agreement with the observations. In section 6 we conclude and give directions for future 
research. 



2 MATHEMATICAL PRELIMINARIES 



In this part we develop formally the spectral decomposition of the elastic wave operator in 
a stratified half-space. We explain how a complete family of generalized eigenvectors of the 
elastic wave equation can be constructed, based on scattering theory. The term "generalized" 
emphasizes the fact that the eigenfunctions are not square integrable, like plane waves. A key 
result is that any square integrable vector field can be expanded over the set of eigenvectors 
of the elastic wave equation, which is an important generalizatio n of the Fourier theorem. 
Detailed calculations and validation of the theory are presented in 




2.1 General properties of the elastic wave operator 



In this part, we summarize results of 



Dermeniian & Guillotl (|1985h and ISecherl (|1998ft . We 



consider the following elastic wave equation in R_, which denotes the semi-infinite medium 
z < 0. The Lame parameters A, fx, and the density p are sufficiently well-behaved functions 
of the position vector x. For instance, we allow these functions to be discontinuous across 
smooth boundaries, but we assume that they are bounded away from zero. The elastodynamic 
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equation reads: 

d tt u(t) + Au(t) = 0, (1) 
where the operator of elasticity acts on the vector wavefunctions as follows: 

1 



Auj(x) 



■dj (X(x)d k u k {x.)5ji + 2//(x)ejj(u)(x)) , 



where ey is the usual strain tensor: 



ey(u)(x) = - (diUj(x) + djUi(x)) . 



(2) 



(3) 



According to Dermenjian and Guillot (1988), the operator A is positive and self-adjoint in 
the Hilbert space with scalar product: 



(u, v) = J (o(x)u* • vd 3 x. 



(4) 



We will not dwell on the precise definition of the domain of A. Physically speaking, the choice 
of the domain ensures that the wavefunctions have finite energy The spectrum of A, i.e. the 
set of generalized eigenvalues of the operator, is the positive real axis and we assume that 
the Lame parameters and the density are sufficiently well-behaved so that the spectrum is 
absolutely continuous (IReed &: Simonlll980l ). In what follows, we briefly explain how to build 
a compl ete set of eigenvector s in stratified media using general results from scattering theory 



in 3-D (Reed fc Simon 



197S 



Ramm 



19861 ). To carry out this task, we first need a basis of 



eigenvectors for the elastic operator in infinite homogeneous space. 



2.2 Eigenvectors of the elastic operator in homogeneous space 

In a homogeneous infinite space the elastic operator takes the well-known form: 

Ao u = -^±^V(V • u) + x V x u. 
Po Po 

One can verify that plane P waves: 

e£(r) 



P" 



ipr 



^ 2 (2vr)3/2 



(5) 



(6) 



are eigenvectors of Ao with generalized eigenvalue p 2 aQ, where ao denotes the P wavespeed, p p 
is a unit vector in the direction of p = (p x ,Py,Pz)- The term generalized is used to emphasize 
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the fact that the eigenfunctions are not square integrable. The notational convention used in 
Equation ([6|) is as follows: the superscript refers to the wave type and the subscript serves 
as a label for the eigenfunction. In the case of a continuous spectrum, the label will usually 
be continuous or sometimes mixed, i.e. discrete+continuous. Similarly, one can introduce two 
independent 5-wave eigenvectors: 

~sh,sv 

sh,8V( \ P is r 

s 1 ' 1/2 /-o \3/2 ' K ' 

where p sh is a unit vector perpendicular to s and the vertical direction z, and p sv is a unit 
vector perpendicular to s and contained in the (z, s) plane. The corresponding eigenvalue is 
/3qS 2 where f3o denotes the shear wave speed. Note that the direction z is for the moment 
arbitrary but will be later taken perpendicular to the stratification. The eigenvectors 
are orthogonal and have the correct continuum normalisation in the sense of the scalar product 
P| defined over M 3 . Let us verify formally the orthogonality of say e^ h and ep: 

/f.sh fj) 
Po(2vrr , . 

R3 (p) 

= p sh -p p S(p-s). 

The delta function requires that p = s in which case the scalar product of the polarization 
vectors vanishes. Thus the scalar product (e^lep 3 ) is identically zero. Let us now show formally 
that our set of eigenvectors is complete, i.e.: 



I<J(r 



Po 



eP(r)eP(r')U 3 p + Jef (r)ef {r')U 3 s + Je s J(r)e s J(r')U 3 s', (9) 



where I is the identity operator in polarization space: Ijj = 6ij and the symbol t means 
hermitian conjugation (conjugate transpose). In equation Q, the normalization prefactor 
Pq 1 is a consequence of the definition of the scalar product (JH). The sums in the right-hand 
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side of equation ([9]) can be written as a single integral: 

e£(r)e£(r')td 3 p + lef (r)ef (r')U 3 s + /e^(rK?(r')W = 



(2vr)^p J v ' 

R 3 

* /e*P(-Vp = (10) 

R 3 

because the set (p p , p sh , p sv ) forms an orthonormal basis in polarization space. We have 
therefore verified the formal completeness relation. 



2.3 Eigenvectors of the elastic wave equation in a stratified half-space 

To construct a complete set of eigenvectors in a stratified Earth, we will apply scattering 
theory. One key idea which can be found e.g. in Reed &: Simon |l979 ) is the following. Let 
us assume that a complete set of eigenfunctions e^ of the wave operator Ao in a reference 
medium is known. Now, perturb the reference medium by introducing a localized scattering 
body to obtain a new operator A. Then it is possible to show that the eigenvectors e^ of 
the perturbed operator are obtained by considering the scattering of e^ by the heterogeneous 
body. More precisely, e^ is the solution of the Lippman-Schwinger equation: 



e fc = e° + G A'e fc , 



where G° is the retarded Green function associated with the operator Aq. The second term 
on the right-hand side of equation (jlip represents the scattered wave which results from 
the interaction between an unperturbed incoming wave and the scatterer represented by A', 
a perturbation of Ao- The set of eigenvectors constructed in this way is orthogonal and 
properly normalized (but not necessarily complet e ) . Sim ilar results have been proved for 



wave scattering by infinite boundaries by 



Ramml (|1986l ). The construction of a c omplete 



Wilcox 



basis of eigenvectors for fluids with arbitrary stratification has been conducted by 
(|1984l ). This author demonstrates how the complete set of eigenfunctions can be constructed 
from th e scattering of plane waves i ncomin g from infini ty. The method of Wilcox has been 



used by 



Dermenjian Guillotl (|1985l ) and ISecherl (|1998l ) to construct rigorously a complete 



basis of eigenfunctions of the elastic operator in a homogeneous half-space with a traction 
free boundary. Applying the same technique, we can exhibit a complete set of eigenvectors 



8 L. Margerin et al. 

for a stratified half-space. Com plete justification of the formulas below are presen ted in a 



Maupin 



separa te paper (|Margerinll2008l ). For an alternative approach the reader is referred to 
(|1996l ). We assume that as z — > — oo, the density and Lame parameters tend to p^, Aqo and 



Hoo, respectively. The decoupling between SH and (P,SV) motions defines two fa milies of 



Wilcox 



soluti ons. The case of SH waves is completely similar to the acoustic case described by 
(1984) and will not be discussed. For the P — SV case we take advantage of the translational 
invariance of the problem and look for solutions of the wave equation which have the following 
form as z — > — oo: 



1 



(2tt)3/V 



m (pV p ' r + r pp pV p ' r + r P s p sv e ls r ) 



where the wavevectors p, p r , and s have the following properties: 



(12) 



P = (Px,Py,Pz), (Px,Py) £K , Pz 
P r = (Px,Py,-Pz), 
S = {px,Py,-^Jp\ (a&JMc ~ 1) +PWo /P'L) > p\ = P^ + Pj- 



(13) 
(14) 
(15) 



As before, p r is a unit vector in the direction of p r and p sv is a unit vector orthogonal to 
s, contained in the (z, s) plane. The P and S wavespeeds at infinity are denoted by and 
Poo, respectively. In equation (|12p . the dependence of the eigenvector on the wavenumber of 
the incident wave is emphasized. This equation can be understood as follows: a generalized 
eigenvector of the operator A is the sum of an eigenvector of Ao , plus the wave back-scattered 
by the layered structure. The reflection coefficients are determined by the continuity of traction 
and displacements across interfaces, and by the vanishing of tractions at the free surface. 
Metho ds to calculate the re fl ection matrix R of the layer stack are described in e.g. 



Kennett 



(j 19831 ) and lAki Richards! (|2002l ) . Equation (|12p defines eigenvectors of A with eigenvalue 



a^p 2 . Another family of solution corresponds to incident SV waves from infinity 



-oo): 



(2vr) 3 / 2 p ( 



1/2 



( P s V s - r + r sp pe ip - r + r ^p^ e isr - r ) 



(16) 



where the wavevector p must be carefully defined, depending on s = (s x , s y ,s z ), s z > 0: 



P 



- 1) + 



>s nv^ 

< S z < Si 



(17) 
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where s?, = s 2 + Sy. The eigenvectors (jX6j) have eigenvalues J 1 = P^s 2 . Note that beyond the 



critical angle for SV — P mode conversion - < s z < sy y — 1— , the reflection coefficients 
as well as the polarization vector p become complex. The reflected P wave is evanescent in the 
lower half-space, and decays exponentially with depth. The complex polarization vector carries 
information on the particle motion of P waves in the half-space. Since the vertical component 
has a — 7r/2 phase shift with respect to the horizontal component, the evanescent P wave is 
in prograde elliptical motion. Equations (jl2H - A16H define the generalized eigenfunctions that 
have oscillatory behavior at infinity. However, they do not form a complete set. Similar to 
bound states that decay rapidly away from a scattering body, the set of eigenvectors must 
be complemented with surface waves, that have their energy localized near the traction- 
free surface. In the case of surface waves, the set of eigenvalues is continuous but organized 
along discrete dispersion branches. The dispersion relation of the surface wave modes can 
be determined by finding the poles of the reflection coe fficients located on the positive real 



axis. 



Numerical procedures are described in the books of 



Kennettl (jl983 ) and 



Aki fc Richards 



(|2Q02| ) . The Rayleigh wave eigenfunctions have the form: 



<,n(r) 



pI + pI 



pI+pI 



(pih™) 



J(pxx+p y y) 



n e N 



(18) 



where d>f n and d>J n denote the horizontal and vertical components of the n Rayleigh 
eigenf unction, and obey the following normalization relation: 



dzp(z) <^ p j(z) 



+ 



1. 



(19) 



These eigenvectors have eigenvalues uj^ = c n (py) pjj, where c n denotes the phase speed of mode 
n. Note that in equation (|18|) . the surface wave eigenvector is labeled simultaneously with the 
continuous indices pn = (p x ,Py) and the discrete mode branch index n. This means that if 
we fix the horizontal wavevector pn and look for solutions of the wave equation that verify 
simultaneously the traction-free condition and the vanishing of traction and displacements at 
depth, this can only occur at a discrete set of eigenfrequencies denoted by u>i, • • • , u n , ■ ■ ■ . Note 
that in applications, we will work the other way around: the central frequency of the signal is 
fixed and one looks for poles of the reflection coefficients related to plane waves with apparent 
velocities smaller than the shear wave velocity in the lower half-space. When a candidate 
mode has been identified, the Rayleigh quotient is used to calculate its group velocity and as 
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an accuracy check. In this formulation, it is implicit that the n th mode appears beyond the 
cut-off frequency uj^. 



3 EQUIPARTITION THEORY 

3.1 Formal definition of equipartition 

Following 



Shapiro et al. 



(200QJ), we introduce the compression and shear deformation energies: 



W p =±p a 2 (V-u) 2 
W s =~Po/3 2 (V x u) 2 



(20) 
(21) 



where ojq, 00 an d Po denote the c qmpressi o nal, shear v eloci ty and the density at the receiver 



and u is the displacement vector. 



Weaver! (1983 , 



19901 s ) and 



Rvzhik et al 



(|1996T ) showed that 



in a heterogeneous infinite elastic medium, the P to S energy ratio stabilizes at large lapse 
time: 

(W s ) _ 2al 

w 



lim 



(22) 



t->+oo (WP) 

where (.) is a statistical averaging over the configuration of scatterers in the medium. This 
result means that in a heterogeneous medium, multiple scattering creates a partition of com- 
pressional and shear energies, which is reached after the waves have encountered a sufficiently 
large number of scatterers. This is the physical interpretation of the limit t — > +oo. Remark- 
ably, the equilibration ratio is completely independent of the nature of the scatterer and of 
the source type. 

The result (|22p gives no information about the time needed to reach equili b rium . The 



significance of the limit t 
experimentally by 



+oo was a ddres sed theoreti c ally 



Malcolm et al. 



elastic radiative transfer equation, 



(2004 ) and 



Pau 



Margerin et al 



et al 



by 



Margerin et al 



(|2000T ) and 



(2005). By solving numerically the 



2000 ) showed that the asymptotic behavior 



occurs after a few mean free times only, where the mean free time denotes the average time 
between two scattering events. These authors also showed that the equilibration time depends 
on the type of scatterers and the type of source. Therefore, the dynamics of the stabilization 
process contains information on the medium heterogeneities. It is also important to realize 
that equilibration is different from equipartition. Equipartition demands that all the modes 
be represented with equal energies whic h in turn implies that the flow of energy i s isotropic 



at lea st away from medium boundaries (jMalcolm et al 



20041 ). In a recent study, 



Paul et al. 



(|2005l ) have shown that equipartition implies equilibration but that the converse statement 
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is wrong. In particular, Equation (|22|) applies even when the energy flux distribution is still 
strongly anisotropic. 



It is mathematically conv enient to define equipartition as a white noise distributed over 



all the modes of the system ([Weaver 



1982; 



Hennino et al. 



2001 



. In order to illustrate this 



definition, let us derive the W s /W p ratio at equipartition in an infinite homogeneous 3-D 
medium using the eigenvectors ©-([T]). Completeness implies that the complex vector field u 
can be expanded as follows: 



u(t,r)= K(r)e— ^d s p + I af ef (r) + af ef(r 



(23) 



where lo t 



(*0\/pI +Pl+P 2 z , Us 



Po\/ s x + s y + s 1- Equation ([23]) defines the analytic 



signal associated with the measured displacements. At equipartition, the amplitudes a p , ... 
are uncorrelated random variables with zero mean and equal variance. In practice, they are 
slowly varying functions of time in the frequency band of interest. In order to investigate the 
energy content at time t in the coda, it is convenient to introduce the Wigner distribution of 
the wavefield: 



E s (t, t, r) = -p /?o (V x u(t + r/2, r) • V x u(t - r/2, r)*) , 



(24) 



where * denotes complex conjugation. In equation (|24p . the brackets denote an ensemble 
average. Assuming that the coefficients a s , a p are governed by a quasi-stationary white- noise 
process: 

(a s (t - r/2)a*,(t + r/2)) = a 2 {t)8(s - s'), (25) 



one obtains: 



E s (t,r,v) = y o /3ia\t) ^ /|V x e™(r)| V^d 3 s. 

„ 1, „„. " 



(26) 



m=sh,sv^ 

In seismological applications, the function a models the slow decay of the energy envelope in 
the coda. The Fourier transform of the Wigner distribution gives the power spectral density 
at frequency u>q: 



E,(t,uo,r) = ^ / L\JI . t. r)c'~" T dr 



(27) 
(28) 



m=sh,sv 
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The assumption of quasi-stationarity demands that the function a does not vary significantly 
during one cycle 2-k/ujq, where ujq is the central frequency of the analyzed signal. Using equa- 
tion ([7]) , one can show that the squared modulus of the rotational part of the eigenvectors of 
a homogeneous full-space is independent of position: 



V x ef(r) 



^ (29) 



Introducing spherical polar coordinates to perform the wave number integral in equation (|28|) . 
one finds the total shear energy spectral density: 



F n ; ,_ ^ 2 (tH ,.-, m 



Using the definition of the compressional energy (|20p and following similar lines, one obtains 
the total compressional energy spectral density: 

£ p (i,WoJ - 77; — r 0—3, (31) 

from which the result (1221) follows. 



3.2 Equipartition in a stratified half-space 



In the case of a stratified half-space, calculations have to be performed numerically and in 
general, the energy ratios depend on depth. Near the boundaries the coupling between P and 
S waves complicates the physical interpretation of t he decomposition of t he field energy into 



Shapiro et al. 



(|2000l ) have shown that, 



shear and compressional components. For instance, 
depending on the incidence angle at the free surface of a homogeneous half-space, a shear 
wave can generate a large amount of compressional energy and vice-versa. In addition, the 
sum of compressional and shear energies does not equal the total elastic deformation energy. 
However, the ratio of shear to compressional energy, as defined above, constitutes a marker of 
the wave content of the seismic coda. We will also consider the vertical to horizontal kinetic 
energy ratio. This quantity can be measured without major difficulties on a rather broad 
frequency range and offers additional information. To illustrate the computational approach, 
let us for instance consider the contribution of the n th Rayleigh mode to the vertical kinetic 
energy E™. Using equation (|18p and the property (|25p . one finds: 

P (z)ula 2 {t) 



8tt 2 



0(pii ,»)(*) 6 { u ° ~ w pii) ^Pb 



(32) 
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where u p ,, = c n (py)p||, and the vector py has components (p x ,Py,0)- Introducing cylindrical 
coordinates, the wavenumber integral can be performed to obtain: 



c n u n 4ir 



(j>u=UJo/Cn,n) 



(33) 



where u n denotes the group velocity of mode n at frequency ujq. The transition from equation 
([32]) to ([33]) uses the following decomposition of the delta function: 



5 (uj - u p A 



1 



CnJ 



(34) 



Equation ([33]) agrees with results of 



Tregoures &: van Tiggelenl (|2002l ) who developed a trans 



port theory for elastic waves in a plate. The interpretation of equation ()33[) is as follows: in 
the seismic coda, the contribution of the n th Rayleigh mode to the vertical kinetic energy is 
proportional to the global density of states of the mode and the local value of the wavefunc- 
tion squared. Thi s product is often referred to as the "local density of states" in the literature 



( Economou 



2003 ) . 



It is also instructive to consider the expression of the kinetic energy for the generalized 
eigenfunctions representing body waves in the lower half-space. For simplicity, we consider the 
kinetic energy of vertical vibrations E% for the generalized eigenfunctions (I12|) corresponding 
to P waves incident from below the sub-surface structure, together with their reflections. To 
facilitate the comparison with the case of Rayleigh waves, we can rewrite the generalized P 
wave eigenvector as follows: 

e ip x x+ip y y 



Y (.p\\>p*n ; 



2tt 



(PlbPz) 



(35) 



i.e. we split the horizontally and vertically propagating parts. The decomposition (j35|) high- 
lights the similar roles played by the surface wave mode index n and the vertical wavenumber 
p z . Using equations (fT2l) and ([25]) . we obtain the following expression: 



8tt 2 



dpz 



^( P| |, Pz )( z )| 5 ( w o-^( P|| , Pz ))dp||, 



(36) 



where ui 



\jVi +Py+Pz, and 

^fpihPz) deno * es * ne vertical component of the gen- 
eralized P-wave eigenfunction, which depends explicitly on the wavenumbers pn and p z . The 
main difference between equations (|36|) and ([32]) is the substitution of the discrete mode 
branch n with the continuous vertical wavenumber p z . After introducing cylindrical coordi- 
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nates, integrating over the vertical wavenumber and the azimuthal angle, one obtains: 



dp\\ 



p\\ 



■P\\ 



'(Pll,p*) 



(37) 



where is the P wave speed in the underlying half-space. Upon making the change of 
variable py 



^-sinS, the last equation simplifies to: 



(t,u) ,z) 



n/2 











2 


1 d6 sin 9 







(38) 



(pn=^ sine '^=^ cose ) 

Evaluated at z = in t he case of a ho mogeneous half-space, equation ([38j) reduces to similar 



expressions obtained by 



Weaver! (|1985l ). based on physical reasoning. Note that away from the 



free surface, the eigenmodes of the system do not reduce to pure P or S waves but consist 
of a superposition of incident and reflected body waves. The result (j38|) can be interpreted 
as follows: the contribution of generalized P eigenfunctions to the kinetic energy of vertical 
vibrations is proportional to the density of states of P waves in the lower half-space and to the 
squared P-wavef unction, averaged over all possible directions of P-waves incident from below 
the structure. Note that the last integral (i38j) also depends on depth and, in general, presents 
complicated oscillations. This is illustrated in Figure [T] where kinetic and deformation energy 
ratios are plotted as a function of depth in a homogeneous 3-D Poisson half-space. The largest 
fluctuations of the deformation and kinetic energy ratios occur in the vicinity of the free surface 
and are mostly caused by the depth dependence of the Rayleigh wave eigenf unction. When 
averaged over one shear wavelength \ s at a depth greater than 3 X s , the energy ratios reduce 
to the expected values for a homo geneous infinit e medi u m. Our calculations are in perfec t 



agreement with previous work by 



Hennino et al 



(|2001I k iTregoures k van Tiggelenl (|2002h . 



based on a locked-mode method. The present approach is generally valid for an arbitrarily 
stratified half-space. Experimentally, the V 2 /H 2 kinetic energy ratio is fairly easy to measure 
over a broad frequency range but can be strongly affected by the local topography. The ratio 
of deformation energies W s /W p is more difficult to evaluate but is independent of the choice 
of local coordinate system. It is therefore expected to be less sensitive to geometrical effects. 
In the next section, we discuss the measurement of the deformation energy ratio on a dense 
array. 
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4 STABILIZATION OF DEFORMATION ENERGY RATIOS IN THE 
SEISMIC CODA 



4.1 Data set and pre-processing 

The data we use were collected by a temporary array that was deployed at the Pinyon Flats 
Observatory (PFO) in 1990. This array operated as part of the Incorporated Research Institu- 
tions for Seismology (IRIS), and was located in Southern California between the San Jacinto 
and the San Andreas fault, as shown on the location maps in Figure [2] 

The array was composed of 58 sensors arranged in a six-sensor by six-sensor grid with 
200m long arms. The sensors were 2Hz L22-D geophones with an acquisition sampling fre- 
quency of 250Hz. They were triggered by two borehole sensors located at 275m and 150m 
depth, respectively. A total of 300 events was recorded, among which 140 were located with 
good accuracy. In our study, we have used 10 events whose source parameters are described 
in Table [TJ These earthquakes have epicentral distances smaller than 50km and magnitudes 
higher than 2. We chose these events because they exhibit a pronounced coda with a high 
signal to noise ratio. We took care of filtering out the frequency components above 40Hz be- 
cause the signal exhibits significant instrumental noise around 60Hz. Because the estimate of 
the diverg ence and curl of the wavefield requires the application of the traction-free boundary 



condition (IShapiro et al 



2000 ). we must first define geometrically the free surface. For this 



purpose we chose the plane which best approximates the location of stations in a least-squares 
sense. The data were rotated in a new reference system whose vertical axis is directed per- 
pendicular to the optimal plane. The signal components were detrended and normalized to 
equal root mean square amplitude at each station i: 



Vt G {1, 



,n} 



u z (ri,t) 2 dt = u 2 z 
u ew (ri,t) 2 dt = u 2 v 
ju ns (ri,t) 2 dt = u 2 ns 



(39) 



where the bar denotes the root mean square. This normalisation procedure corrects for possible 
instrumental response differences and yields more reliable estimates for the field derivative. 



4.2 Method of estimation for P and S energies 

To estimate the deformation energies W p and W s , the displacement gradients diUj of the 
wavefield need to be measured. We will first examine the problem of estimating the horizontal 
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derivatives. These will be subsequently used to obtain the vertical derivatives. Using three or 
four stations as shown in Figure [3l the derivative of the wavefield in two linearly independent 
horizontal directions can be estimated, which suffices to recover the horizontal components of 
the gradient vector by application of the Taylor expansion formula: 



U( r 2) 



n(ri) + (r 2 - ri).Vu + O ( |r 2 — ri 



(40) 



In practice, the estimate of the spatial derivative requires a great coherence of the displace- 
ments at two near-by stations. Applying a finite difference approximation to equation (|40p , 
one finds 

duj ^ Uj(ri + dj) - Uj(ri) 
dxj ~ 



d 



(41) 



where d is the inter-station distance. Ui(r\ + d.j) and Ui(ri) correspond to the measured 
displacements at positions r 2 = ri + dj and rj, respectively. The approximation (I41|) will 
be correct provided d/X < C 1 a nd its application to the data set will be critically examined. 



Following 



Shapiro et al. 



(2000), we use the free surface boundary condition to express the 



derivative in the z direction in terms of horizontal derivatives: 



d z u x = 


-d x u z 


8 z Uy = 


-d y u z 


d z u z = 


-(1- 



(1 - 2f3 2 /a I ) [d x u x + d y u 

It follows that the divergence and curl of the wavefield can be expressed in terms of derivatives 
along horizontal coordinates exclusively 



a 

(2d y u 



V • U (9xU X + OyUyJ . 

V x u 



y u z ) 



2,d x u z , d x Uy 



d y u x ) 



(42) 
(43) 



We now study the validity of our estimates of spatial derivatives in the data. Practically, 
we approximate d x u x with the finite-difference formula (I41j) . The estimate is based on the 
assumption that the two signals are very similar, i.e., the change of waveform between two 
stations is infinitesimal. Using the converse assumption we propose a method to determine the 
limit of validity of the finite difference approximation. When two signals are poorly correlated 
the following result applies: (itj(l)iij(2)) <C uf. Therefore we can obtain the following estimate 
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of the derivative for incoherent signals 



<(u,(2) - Ui (l))2) _ ^K 2 (2)> + («?(!)) - 2(^(1)^(2)) 




(44) 



Therefore we expect to find al/d behavior when the signals at 2 stations are poorly correlated, 



inter-station distance d. In the 5 — 7Hz frequency band, and for distances shorter than 50m, 
the evaluation is unstable because the time shifts between the two waveforms becomes nearly 
equal to the sampling period. For distances greater than 150m, we roughly find the predicted 
1/d behavior (see equation (|44p ). For inter-station distances ranging from 50m to 150m, the 
finite difference curve shows a plateau, indicating that the measurement is stable and reliable. 
In order to avoid possible biasing by extreme values that sometimes appear, we define the 
experimental derivative as the median of all available finite difference values. 

4.3 Stabilization of shear to compressional energy ratios 

Compressional (W p ) and shear (W s ) energies have been measured for each event and for 
different frequency bands. Figure [5] shows a typical example of the band-pass filtered signal 
between 5 and 7 Hz, the decay of energies with time, and the dynamics of the ratio between 
shear and compressional energies. For this event a stabilization of the energy ratio is clearly 
visible in the coda. While the total energy itself decreases by two orders of magnitude, the 
ratio shows remarkably weak fluctuations. In addition, we infer that the stabilization occurs 
very rapidly, only a few cycles after the arrival of the direct S wave. When the compressional 
energy W p reaches the background noise level, the energy ratio again fluctuates randomly 
and stabilization disappears. This observation indicates that the wave contents of noise and 
coda are rather different. In a sense, the coda is easier to interpret because it is likely to 
contain an equipartition mixture of all the modes of the system. Because fluctuations are 
weak, we define the stabilization ratio as the average of W s /W p over a time window starting 
one second after the direct S arrival and ending when W p is of the order of the noise level. 
The equilibration has been observed for the ten earthquakes provided that frequencies lower 
than 3 Hz and greater than 9 Hz have been removed. The low frequency coda is difficult to 
measure because of the poor sensitivity of the sensors. At too high frequencies the diffuse 
field is not sufficiently coherent to obtain reliable estimates of the derivatives. The results are 
shown in Figure El where we illustrate the fact that the equipartition ratio is independent 



i.e. sufficiently far apart. Figure d] shows the finite difference 



Ui(ri+dj)—Ui(ri) 
d 



clS cl function of 
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of the source parameters. This favors an interpretation of the observed stabilization as a 
consequence of multiple scattering. However, the equilibration ratio equals 2.8 +/- 0.4, which 
is much lower than the ratio 7.2 predicted at the surface of a homogeneous Poisson half- 
space (Figure [1]). A likely explanation is the non-uniformity of vel ocities in the upper c rust. 



The seismic propert ies under the PFO array have been studied by 



Vernon et al. 



Fletcher et al. 



(|1990h and 



(1998). These authors have shown that the first 50 to 70 meters at Pinyon Flats 
are composed of weathered granite with low seismic velocities. At greater depth, intact granite 
with high wavespeeds is found. Therefore, we speculate that the anomalously low energy ratio 
at PFO is caused by the presence of low-velocity layers. We further explore this hypothesis 
by applying equipartition theory in stratified media. 



5 MODELING ENERGY RATIOS IN THE CODA 

he s tructure under the PFO array was explored by 



Fletcher et al 



(|199d ) and 



Vernon et al. 



( 19981 ) . We first study a simplified model of the velocity and density profiles, composed of 



a thin layer (65 m thickness) with relatively low velocities (2.77 km/s for P-wave and 1.6 
km/s for S'-wave) overlying a half-space with typical crustal velocities (5.2 km/s for P-wave 
and 3 km/s for S'-wave). The elastic medium is poissonian everywhere and the density is 
uniform. The properties of the different seismic models studied in this paper are summarized 
in Table [2j In Figure [7J the frequency dependent theoretical ratio between shear (W s ) and 
compressional (W p ) energies at equipartition is represented. At low frequency, the waves are 
insensitive to the thin low-velocity layer and the 7.2 ratio of a homogeneous Poisson half-space 
is recovered. At low frequency, the fundamental Love wave can be completely neglected, and 
the largest contribution to the S and P deformation energies is made first by the Rayleigh 
waves, and second by the generalized eigenfunctions representing body waves in the lower 
half-space. At very high frequency the W s /W p ratio again converges to the half-space value. 
This is not so surprising, since we may expect the diffuse wavefield to be rather insensitive to 
the deep velocity contrast. The pronounced drop of the energy ratio to values smaller than 
4 around the fundamental resonance frequency of the layer is largely due to the change of 
the Rayleigh wave ellipticity and the increasing importance of the fundamental Love mode 
near the resonance frequency of the layer. At high frequencies, the generalized eigenfunctions 
corresponding to body waves in the underlying half-space play little role in the equipartition 
ratio. The high-frequency oscillations are mostly due to the interplay between Rayleigh and 
Love waves. 

Because the deformation energies are accessible in a narrow frequency band only, we have 
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also evaluated the ratio of vertical (V 2 ) to horizontal (H 2 ) kinetic energies in the 5-25 Hz 
band. The measurement of this ratio is straightforward and offers additional information on 
the local structure. The experimental results are plotted in Figure [5] and show a clear fre- 
quency dependence. The ratio decreases from 0.5 at frequencies around 4 Hz to about 0.1-0.2 
in the 6-12Hz frequency band, then rapidly increases to 0.8, and oscillates around this value 
in the 15-25Hz frequency band. While a simplified 1-layer model (model 1 in Table [2]) gives a 
clear qualitative explanation of the observations at PFO, it is not quantitatively satisfactory. 
By a process of trial and error and using a-priori information provided by previous studies, 
we tried to find 2-layer models that better match the observation. The velocity and density 
profile of a typical 2-layer model is given in Table [2 The velocities in the bedrock and i n 

f|1990l) . 



Fletcher et al 



the deepest part of the weathered zone are identical to those given by 
According to these authors the near-surface seismic properties are more variable and less well 
constrained. Therefore we explored various thicknesses and velocities for the top layer. Re- 
garding the density, little information could be found in the literature. We adopted a typical 
value of 2.7 for the intact bedrock and allowed the density to decrease to 2.2 in the weathered, 
less consolidated layers. Observed and theoretical energy ratios for a 2-layer model that fits 
reasonably well the W s /W p and V 2 /H 2 ratios up to 15 Hz are shown in Figure El In the sim- 
ple 2-layer model, the small V 2 /H 2 ratio at low frequency is due simultaneously to the nearly 
horizontal polarization of the fundamental mode Rayleigh wave and the increasing role played 
by the fundamental Love wave. For a thorough study of the freq uency-dependent ellipticity of 



the R ayleigh wave in layered structures, the reader is referred to iMalischewskv fc Scherbaum 
(200J). In model 2 (see Table [2]), at frequencies higher than 5 Hz, the surface waves trapped 
in the low-velocity layers largely dominate over body waves coming from below the structure. 
The sharp increase of the V 2 /H 2 ratio at high-frequencies is again explained by an interplay 
between Love and Rayleigh waves. Around 15 Hz, the fundamental mode Rayleigh wave be- 
comes strongly vertically polarized and the weight of Love waves in the equipartition ratio 
suddenly drops. The two effects combined explain the peak value observed in the data. Cal- 
culations based on Rayleigh waves only, predict fluctuations of the V 2 /H 2 ratio much bigger 
than observed. This demonstrates that the contribution of Love waves is of fundamental im- 
portance. At frequencies higher than 15Hz, the measured kinetic energy ratio oscillates and 
departs significantly from the model calculations. This may indicate a departure from the 
simple layered structure. 



By exploring a large number of models, we found that the position of the steep rise from 
the low to high frequency behavior is extremely sensitive to the the thickness and velocity 
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of the upper layer. This is illustrated in Figure [9] where the frequency dependent vertical to 
horizontal kinetic ratio has been measured at one station located at the end of one of the 
two orthogonal arms (see Figure [2]). Figure [9] also shows the calculations of the frequency 
dependent ratio between the vertical and horizontal kinetic energies for model 3, Table [2j 
This model is completely similar to model 2, except in the first 15 meters where the velocities 



are slightly different from those proposed by 



Fletcher et al. 



(|l990l ). This model fits reasonably 



well the observations over the whole frequency band. The comparison of Figures [8] and [9] 
demonstrates that the frequency dependence of the ratio between vertical and horizontal 
kinetic energies measured in the coda contains information on the local velocity structure. In 
the case of PFO, our study supports the idea that the velocity profile presents simultaneously 
velocity gradients at depth and thin very low-velocity layers at the surface. 



6 CONCLUSION AND OUTLOOK 



We have observed the stabilization of ratios of deformation and kinetic energies in the coda of 
small earthquakes, at PFO. The stabilization phenomenon is interpreted as a sign of multiple 
scattering of waves in the crust at high frequencies. The ratio between the shear and compres- 
sional energies can be accurately estimated in the 5-7 Hz frequency ba nd and was observed 
to be much lower than the one found in Mexico bv lShapiro et al.l (|2000l ). To understand this 
observation, we have developed a theory of equipartition in arbitrary layered elastic media 
based on the spectral decomposition of the elastodynamic operator. We have shown that the 
decrease of the S to P deform ation energy ratio may be explained by the low-velocity sub- 



surface layers. Like in Mexico ([Shapiro et al 



2000; 



Hennino et al 



2001 



), the stabilization of 

various energy ratios occurs shortly a fter the S wave arriva l. Such a rapid stabilization has also 
been found in numerical simulations (IMargerin et al.ll2000l ). typically after a moderate number 
(3-4) of scattering events. This suggests that the mean free time -the typical time between 
two scattering events- is very small, of the order of a few seconds only, and supports the idea 
of a highly heterogeneous crust in California. Our results on the frequency dependence of the 
ratio of vertical to horizontal kinetic energies sound reminiscent of the so-called Nakamura's 
technique which relie s on the fr equency dependence of the H/V spectral ratio of ambient 



noise vibrations (see lBardlll998l . for a review). This technique usually makes use of the low- 
frequency H/V peak to retrieve the resonance frequency of sedimenta ry deposits. The physical 



Bonnefoy-Claudet et al.1 



bases and limitations of the H/V method are still actively debated. 
(2006) present a detailed study of a simple 1-layer configuration which illustrates the com- 
plex interpretation of H/V measurements depending on the spatial distribution of the noise 
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sources. Because coda waves are composed of an equipartitioned mixture of all modes, they 
are in a sense easier to model than noise signals. We have shown that it is possible to model 
reasonably well the frequency dependent ratio of vertical to horizontal kinetic energies in the 
coda at PFO. In addition to a low-frequency global resonance of the 1-D structure, we find 
that some information on the uppermost layers can also be obtained. Another interesting 
aspect of the coda is the possibility to check the wave content of the signal through the mea- 
surements of several energy ratios that can be compared to theoretical predictions. Further 
investigations of site effect assessments with coda waves are required before reaching more 
definitive conclusions. 
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Origin time 


Latitude 


Longitude 


Magnitude 


Depth 


Source-array distance(km) 


117.12.52.42.392 


34.0487 


-116.3908 


2.3 


-2.2 


48.89 


118.08.32.38.165 


33.8647 


-116.1659 


2.4 


5.1 


39.04 


122.11.35.03.691 


33.4927 


-116.4624 


2.1 


7.0 


13.20 


125.08.10.15.667 


33.5074 


-116.4672 


2.0 


8.3 


11.64 


127.12.40.57.581 


33.8693 


-116.1555 


2.0 


3.7 


40.08 


130.07.23.35.987 


33.6496 


-116.7291 


2.1 


15.7 


25.46 


130.14.25.10.270 


33.1944 


-116.3622 


2.8 


14.5 


47.13 


132.23.54.49.190 


33.9839 


-116.3103 


2.5 


4.8 


43.50 


134.05.05.27.123 


33.5384 


-116.6073 


2.6 


11.8 


16.01 


136.01.14.21.003 


33.4376 


-116.4491 


2.4 


11.1 


19.30 



Table 1. Location and magnitude of the 10 events used in this study. 
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Model 1 


Model 2 


Model 3 






hi = 4m 






ai = 300m/s 






Pi = 15Um/s 




hi = 11 m 


pi = 2200kg/m 3 




ai = 720m/s 


hi = 65m 


fix = 400m/s 


/i2 = 11m 


ai = 2770m/s 


pi = 2200kg/m 3 


a 2 = 900m/s 


Pi = 1600m/s 




/3 2 = 500m/s 


pi = 2700kg/m 3 


h 2 = 41 m 


p 2 = 2200kg/m 3 




«2 = 3100m/s 






P 2 = 1600m/s 


/13 = 50m 




p 2 = 2500kg/m 3 


a 3 = 3100m/s 






03 = 1600 m/s 






p-i = 2700kg/m 3 


aoo = 5200m/s 


ctoo = 5400m/s 


aoo = 5400m/s 


/?oo = 3000m/s 


/3oo = 3000m/s 


/?oo = 3000m/s 


Poo = 2700 kg/m 3 


Pec = 2700kg/m 3 


Poo = 2700kg/m 3 



Table 2. Seismological models of the sub-surface at Pinyon Flats Observatory us ed in this study. The 
different depths and velocities are inspired by the results of IFletcher et al.l (| 19901 ). 
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1 2 3 4 5 

Depth (A s ) 




1 2 3 4 5 

Depth (X s ) 



Figure 1. Depth dependence of energy ratios at equipartition in a homogeneous elastic half-space 
(Poisson solid). The depth unit is the shear wavelength X s . Top: vertical to horizontal kinetic energy 
ratio. Bottom: shear to compressional deformation energy ratio. Note the persistent oscillations at 
depth, which originate from the interference between incident and reflected waves at the free surface. 
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117 116 115 114 113 



Figure 2. Location map of the 1990 Pinyon Flats High Frequency Array Experiment. 58 3-components 
sensors were deployed during 3 months. The upper-left inset shows the location of the largest magni- 
tude local events recorded by the array. The lower-left inset shows the acquisition geometry: a dense 
square of 6 x 6 seismometers separated by 7 meters, and two perpendicular l ong arms composed of 11 
seismometers separated by 21 meters. Reproduced from lVernon et a l. (1998) 



Figure 3. Estimation of the divergence and curl of the wavefield requires the measurement of the 
derivative along two linearly independent directions. This figure presents two valid configurations with 
3 or 4 stations 
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Figure 4. Measurement of the spatial derivative of the wavefield on a dense array of seismometers. 
Dots: estimates of the derivative du z /dx for one event at lapse time t = 32.4s. The mean value of the 
field derivative is indicated by a solid line. The horizontal dotted lines indicate the plus or minus one 
standard deviation range. The + symbols show the estimate of the derivative for incoherent fields. The 
plot illustrates the stability of the estimate of the derivative for stations located at least 50m and at 
most 150m apart. Beyond this distance the field at the two stations become uncorrelated as shown by 
the coincidence between the incoherent and coherent estimates of the wavefield derivative. 
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Figure 5. Observation of the stabilization of deformation energies in the seismic coda. Top : example 
of vertical displacement for a shallow magnitude 2.0 local event. Middle : compressional and shear 
deformation energies as a function of time. Note the logarithmic scale on the vertical axis. Bottom : 
time dependence of the shear to compressional energy ratio. The ratio shows small random fluctuations 
while the total energy decays by two orders of magnitude. Note the large fluctuations of energy ratios 
in the noise following the coda. The horizontal line indicates the time window used to estimate the 
stabilization value. 
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Figure 6. Estimate of shear to compressional energy ratios for 10 local events in the 5 — 7Hz frequency 
band. The stabilization ratio varies little from one event to the other. The shaded region indicates plus 
or minus one standard deviation around the mean value. This feature agrees with the theoretical 
prediction that the energy ratio in the coda is independent from the source parameters. 
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Figure 7. Frequency dependence of the equipartition ratio between shear iW s ) and compressional 
(W p ) energies in a one-layer model of Pinyon Flats Observatory. At a depth of 65m, the shear velocity 
increases by a factor almost 2 (model 1 in Tabled]). At low and high-frequencies, the equipartition 
ratio of a Poisson half-space is recovered. Close to the resonance frequency of the low- velocity layer, a 
significant drop of the ratio W s /W p is observed. 
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Figure 8. Modeling of the measured stabilization ratios with a simple 2-layer model of PFO as 
described in Table [2] (model 2). Top: Observed (shaded region) and modeled (dashed line) shear to 
compressional energy ratio. Bottom: Observed (shaded region) and modeled (dashed line) vertical to 
horizontal kinetic energy ratio. The black solid lines delimit the ±1 standard deviation region around 
the mean value. Agreement is reasonable up to 15 Hz. 
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Figure 9. Modeling of the measured stabilization ratio between vertical and horizontal kinetic energies 
with a 3-layer model of PFO as described in Table [2] (model 3). The shaded area delimits the ±1 
standard deviation region for the data. The dashed line shows the model calculations. 



